Comparing models for background adjustment

The point of this document is to identify an appropriate function to use to adjust for the gradient seen in the EM videos. The two options considered so far are

• A quadratic relationship between the \(x\)-coordinate and the pixed intensity (\(z\)), i.e., \(z=\mu_g+\beta_1x+\beta_2x^2\)
• A nonparametric relationship between the \(x\)-coordiante and the pixel intensity, i.e, \(z=\mu_g+\beta_1x+\beta_2x^2+\lambda\sum_{j=1}^Jf(x-k_j)\) where \(k_j\) are user specified knots and \(f(\cdot)\) is some function such as truncated polynomial of a chosen degree and \(\lambda\) is a penalty parameter

The problem with these methods is that the former is not flexible enough and the latter is a bit too flexible. In particular, the quadratic relationship cannot account for differences in the gradient across the node and pure background. The spline can account for these different gradients, but it would be hard to generalize because there are so many knobs associated with splines (number and location of knots, penalty).

To illustrate, consider frame 4292 of the first video. The top image is the full picture while the bottom two are the pixel intensities as a fucntion of \(x\) and \(y\) for the 5% sampled images.

g1 <- qplot(x,y,colour=Intensity,data=m1_dat)+theme_em()+scale_colour_continuous(low="black",high="white")
g2 <- qplot(x,Intensity,data=m1_10)
g3 <- qplot(y,Intensity,data=m1_10)
grid.arrange(g1,arrangeGrob(g2,g3,nrow=1),nrow=2)

Because the slope of the gradient seems to differ for pixels on the node and pixels in the background, we may consider adding an interaction term with \(y\). In particular we can extend the quadratic model in \(x\) to \[z=\mu_g+\beta_1x+\beta_2x^2+\beta_3yx\;\text{ or }\;z=\mu_g+\beta_1x+\beta_2x^2+\beta_3yx+\beta_4xy^2.\] To see how these functions look for these data, I fit the basic quadratic model and each of the interaction models to the pixel intensity values for the frame above to get a sense for how the gradient adjustment will look.

#Include xy^2 interaction
x2ylm <- lm(Intensity~I(stdize(y))+I(stdize(x*y))+I(stdize(x*x*y))+I(stdize(x))+I(stdize(x^2)),data=m1_10)
x2ysurf_data <- data.frame(x=x2ylm$model[,5], y=x2ylm$model[,2], z=x2ylm$fitted.values)
g1 <- qplot(x,y,data=x2ysurf_data,colour=z,size=I(4))+theme_em()+
  ggtitle(expression(beta[0]+beta[1]~x+beta[2]~x^2+beta[3]~xy+beta[4]~x^2~y))+scale_colour_continuous(low="black",high="white")
#Only xy interaction
xylm <- lm(Intensity~I(stdize(y))+I(stdize(x*y))+I(stdize(x))+I(stdize(x^2)),data=m1_10)
xysurf_data <- data.frame(x=xylm$model[,4], y=xylm$model[,2], z=xylm$fitted.values)
g2 <- qplot(x,y,data=xysurf_data,colour=z,size=I(4))+theme_em()+
  ggtitle(expression(beta[0]+beta[1]~x+beta[2]~x^2+beta[3]~xy))+scale_colour_continuous(low="black",high="white")
#No interaction
xonlylm <- lm(Intensity~I(stdize(x))+I(stdize(x^2)),data=m1_10)
xlm <- data.frame(x=xonlylm$model[,2],y=xylm$model[,2],z=xonlylm$fitted.values)
g3 <- qplot(x,y,data=xlm,colour=z,size=I(4))+theme_em()+
  ggtitle(expression(beta[0]+beta[1]~x+beta[2]~x^2))+scale_colour_continuous(low="black",high="white")
grid.arrange(g3,g2,g1,nrow=1)

Now implement those gradient models with the GMM to see how the predicted pixel intensities compare to the observed values. It’s probably worth checking that I’m computing the predicted intensities correctly, so the predicted intensity \(\hat z\) for pixel \(i\) is \[\hat z_i=\mathbf{x}_i^\top\hat{\mathbf{\beta}}+\sum_{g=1}^G\hat{w}_{ig}\hat{\mu}_g\] where \(\mathbf{x}\) is the appropriate vector of \(x\)- and/or \(y\)-coordinates, \(w_{ig}\) is the probability that pixel \(i\) belongs to group \(g\) and \(\mu_g\) is the group \(g\) mean.

#No interaction
std_quad_reg <- list(p=rep(1/4,4),beta=unname(xonlylm$coefficients[-1]),delta=c(50,100,150,200),sigma=rep(6,4))
std_quad_reg_res <- classify_growth_LR(z=m1_10$Intensity, x=m1_10$x, y=m1_10$y, inits=std_quad_reg, eps=1e-2, maxit=2000,interaction = FALSE)
m1_10$Std_reg <- gmm_zhat(em_out = std_quad_reg_res$em_out,xmat = std_quad_reg_res$X)
std_x <- qplot(x,Intensity,data=m1_10)+geom_point(aes(x,Std_reg),colour=2,shape=1, alpha=.5)
std_y <- qplot(y,Intensity,data=m1_10)+geom_point(aes(y,Std_reg),colour=2,shape=1, alpha=.5)
#xy interaction
init_reg1 <- list(p=rep(1/4,4),beta=unname(xylm$coefficients[-1]),delta=c(50,100,150,200),sigma=rep(6,4))
int_quad_reg_res <- classify_growth_LR(z=m1_10$Intensity, x=m1_10$x, y=m1_10$y, inits=init_reg1, eps=1e-2, maxit=2000,interaction = TRUE, interaction_degree = 1)
m1_10$Deg1 <- gmm_zhat(em_out = int_quad_reg_res$em_out,xmat = int_quad_reg_res$X)
deg1_x <- qplot(x,Intensity,data=m1_10)+geom_point(aes(x,Deg1),colour=2,shape=1, alpha=.5)
deg1_y <- qplot(y,Intensity,data=m1_10)+geom_point(aes(y,Deg1),colour=2,shape=1, alpha=.5)
#xy, xy^2 interaction
init_reg2 <- list(p=rep(1/4,4),beta=unname(x2ylm$coefficients[-1]),delta=c(50,100,150,200),sigma=rep(6,4))
LR_int_res <- classify_growth_LR(z=m1_10$Intensity, x=m1_10$x, y=m1_10$y, inits=init_reg2, eps=1e-2, maxit=2000,interaction = TRUE,interaction_degree = 2)
m1_10$Deg2 <- gmm_zhat(em_out = LR_int_res$em_out,xmat = LR_int_res$X)
deg2_x <- qplot(x,Intensity,data=m1_10)+geom_point(aes(x,Deg2),colour=2,shape=1, alpha=.5)
deg2_y <- qplot(y,Intensity,data=m1_10)+geom_point(aes(y,Deg2),colour=2,shape=1, alpha=.5)
grid.arrange(arrangeGrob(std_x,std_y,top=grid::textGrob(label=expression(hat(mu)+hat(beta)[1]~x+hat(beta)[2]~x^2))),
             arrangeGrob(deg1_x,deg1_y,top=grid::textGrob(label=expression(hat(mu)+hat(beta)[1]~x+hat(beta)[2]~x^2+hat(beta)[3]~xy))),
             arrangeGrob(deg2_x,deg2_y,top=grid::textGrob(label=expression(hat(mu)+hat(beta)[1]~x+hat(beta)[2]~x^2+hat(beta)[3]~xy+hat(beta)[3]~xy^2))),nrow=1)

To see how well these predicted intensity values are translated into growth/non-growth. The pictures below have not be spatially corrected.

#No interaction
m1_10$std_growth <- std_quad_reg_res$Growth
std_grow <- qplot(x,y,colour=Intensity,data=m1_dat)+geom_point(aes(x,y),colour=2,shape=1, data=filter(m1_10,std_growth))+
  theme_em()+scale_colour_continuous(low="black",high="white")
#xy interaction
m1_10$int_growth <- int_quad_reg_res$Growth
int_grow <- qplot(x,y,colour=Intensity,data=m1_dat)+geom_point(aes(x,y),colour=2,shape=1, data=filter(m1_10,int_growth))+
  theme_em()+scale_colour_continuous(low="black",high="white")
#xy, xy^2 interaction
m1_10$xyint_growth <- LR_int_res$Growth
int2_grow <- qplot(x,y,colour=Intensity,data=m1_dat)+geom_point(aes(x,y),colour=2,shape=1, data=filter(m1_10,xyint_growth))+
  theme_em()+scale_colour_continuous(low="black",high="white")
grid.arrange(arrangeGrob(std_grow,top=grid::textGrob(label=expression(beta[1]~x+beta[2]~x^2))),
             arrangeGrob(int_grow,top=grid::textGrob(label=expression(beta[1]~x+beta[2]~x^2+beta[3]~xy))),
             arrangeGrob(int2_grow,top=grid::textGrob(label=expression(beta[1]~x+beta[2]~x^2+beta[3]~xy+beta[3]~xy^2))),nrow=1)

Based on the plots of fitted and observed pixel values as a function of \(x\), it looks like just the \(xy\) variable is needed and not \(xy^2\). So apply the the correction with and without an \(xy\) interaction term to see how they perform on the full video 1 and 3. These videos have been spatially corrected with a radius of 10 and \(\gamma\) value of 0.2.